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Abstract. Finite volume schemes are commonly used to construct approximate solutions to 
conservation laws. In this study we extend the framework of the finite volume methods to 
dispersive water wave models, in particular to Boussinesq type systems. We focus mainly on 
the application of the method to bidirectional nonlinear, dispersive wave propagation in one 
space dimension. Special emphasis is given to important nonlinear phenomena such as solitary 
waves interactions. 



1. Introduction 

The simulation of water waves in realistic and complex environments is a very challenging 
problem. Most of the applications arise from the areas of coastal and naval engineering, but also 
from natural hazards assessment. In this work we will study numerically bidirectional water wave 
models. Specifically, we consider the following family of Boussinesq type systems of water wave 
theory, introduced in [4], written in nondimensional, unsealed variables 



(1.1) 



■rit+Ux + {'nu)x + auxxx - brjxxt = 0, 

Ut + Vx + UUx + CT]xxx - dUxxt = 0, 

where a, 6, c, d G R, 77 = 77(0;, t), u — u{x, t) are real functions defined for a; e R and t > 0. 

a = ^(^^ - l)^, b^\i&'- - c = 1(1 - e'),,, d^^ii- e^)ii - m), 

where <6 <1 and ^, G R. 

Finite volume method is well known for its accuracy, efficiency and robustness for approximat- 
ing solutions to conservation laws and in particular to nonlinear shallow water equations. The 
aforementioned bidirectional models (1.1) are rewritten in a conservative form and discretization 
by the finite volume method follows. Three different numerical fluxes are employed 

• a simple average flux (m-scheme), 

• a central flux, (KT-scheme) [16, 14], as a representative of central schemes, 

• a characteristic flux (CF-scheme), as a representative of the linearized Riemann solvers, 
[10]. 

along with TVD, UNO and WENO reconstruction techniques, [18, 12, 15]. Time discretization 
is based on Runge-Kutta (RK) methods which preserve the total variation diminishing (TVD) 
property of the finite volume scheme, [17]. We use explicit RK methods since we work with BBM 
type systems (1.1) and not with KdV equation which is well known to be notoriously stiff. 

The present text is organized as follows. In Section 1 we present the mathematical model 
under consideration and the context of this study. Then, Section 2 contains a brief description 
of various numerical schemes we implemented. Accuracy tests and several numerical results on 
head-on collision of solitary waves are presented in Section 3. Finally, some conclusions of this 
study are outlined in Section 4. 
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2. Numerical schemes 

In the present section we generalize the finite volume method to systems (1.1) of dispersive 
PDEs. Boussinesq system (1.1) can be rewritten in a conservative like form as follows: 

(I-D)vt + [F(v)], + [G(v)],=0, (2.1) 

where v = (?7, u)^, F(v) = ((1 + r])u,r] + ^u^)'^, G(v) = {auxx,cr]^x), and D = diag (6 9^, d5^). 
The simplest discretization is based on the average fluxes J^™ for F and C/™ for G. For the other 
two choices of the numerical flux T the evaluation of Jacobian is needed. Let A denotes the 
Jacobian of F, then 

( ^ 1 + ^^ 
1 u 



with eigenvalues \i — u± vT+Ty, i — 1,2. It is readily seen, since F is a hyperbolic flux, that 
A can be decomposed as A = LAR thus for the characteristic flux J^'^'^ we have with ii — , 
Si = sign(Ai), 1 = 1,2 

A{W, V) -- 



For evaluating the numerical fluxes G simple cell averages or higher order approximations such 
as UN02 or WENO can be used. For more details we refer to our original research article [9]. 

Remark 1. The discretization of the elliptic operator D is based on the standard centered dif- 
ference. This is a second order accurate approximation and it is compatible with the TVD2 and 
UN02 reconstructions. For higher order interpolation we need to modify the elliptic and flux dis- 
cretization to match the reconstruction's order of approximation. Indeed, the finite volume scheme 
is modified as 

d 



dt 



V,_i + lOV, + V,+i ,^ ^. V,+i - 2V, + V,; 



-{b,d)- 



12 ' ' ' Aa;2 



12 ^ 



where Hi = ■^(-^i+i ^ -^i-J-) + -^{Gi+i ^ Gi^i), is a fourth order accurate approximation. 

Remark 2. In the sequel for the discretization of the dispersive term G we use mainly the average 
numerical flux G"^ defined as G™^i = (a, c) ^ where = "^'^^ A^a"*""^'"^ ■ ^'^ '^'^^^ of higher 

order WENO reconstructions we use the average numerical flux based on the reconstructed values 

of Y," i.e. the flux C = (a,c) ^ —, where Y , , and Y , are reconstructed values of 

Y,. 

2.0. 1. Boundary conditions. In the case of Bona-Smith type systems with flat bottom we consider 
herein only the initial-periodic boundary value problem which is known to be well-posed [1]. 

3. Numerical results 

For the Boussinesq system (1.1) we present first results demonstrating the accuracy of the finite 
volume scheme. Then, we study interaction of solitary waves. 

3.1. Accuracy test, validation. We consider the initial value problem with periodic boundary 
conditions for the Bona-Smith systems with known solitary wave solutions [7] to study the accuracy 
of the finite volume method: 

viO = ?7osech^(AO, 

with 

_ 9 9^-7/9 _ 4(e^-2/3) 

VO - 2 ■ 1-02 , Cs - ^2(i _g2)(g2_ i7^' 

^ _ 1 / 3(fl2-7/9) „ _ / 2(l-92) 

^ — 2 Y (e^-i/3)(e2-2/3) ' ^ — 
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Table 1. Rates of convergence. 
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Figure 1 . Preservation of the solitary wave amplitude and conservation of the 
invariant : G™ flux with Minmod hmiter 



We fix 9^ — 8/10 in the system and an analytic solitary wave of amplitude rjo — 1/2 is used as 
the exact solution in [—50, 50] computed up to T = 100. The error is measured with respect to 
discrete and L°° norms, namely we use: 

Elik) = \\U%/\\U'h, \\U% = (^Ax\U^\^ 
E^{k) = ||[/"|U,oo/||(7°||„,oo, ||t/'||/..oo = max|C/f I, 

i 

where — denotes the solution of the fully-discrete scheme at the time — k At. The 

expected theoretical order of convergence was confirmed for all finite volume methods we presented 
above. Two indicative cases are reported in Table 1 for the average flux and TVD2 implementation 
with MinMod limiter. 

We also check the preservation of the invariant Ii{t) = /g(?7^(a;,t) + (1 + r]{x,t))u'^{x,t) — 
crj^ixjt) — au1{x,t)) dx by computing its discrete counterpart: 

as well as the discrete mass Iq = Ax Vi- Figure 1 shows the evolution of the amplitude and the 
invariant of the solitary wave up to T = 200. The comparison of various methods is performed. 
We observe that the UN02 reconstruction is more accurate while KT and the CF schemes show 
comparable performance. We note that the invariant Iq — 1.932183566158 conserved the digits 
shown for all numerical schemes. In this experiment we took Aa; = 0.1, At = Ax/2. 
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(a) Bona-Smith 
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Table 2. Preservation of the invariant /{' 



3.2. Head-on collisions. The head-on colhsion of two counter-propagating solitary waves is char- 
acterized by the change of the shape along with a small phase-shift of the waves as a consequence 
of the nonlinearity and dispersion. These effects have been studied extensively before by numer- 
ical means using high order numerical methods such as finite differences, [3], spectral and finite 
element methods [2] and experimentally in [8]. In Figure 2 we present the numerical solutions of 
the BBM-BBM system and the Bona-Smith system with 6*^ = 9/11 (in dimensional and unsealed 
variables) along with the experimental data from [8] . The spatial variable is expressed in centime- 
ters while the time in seconds. The solutions were obtained using the CF-scheme with UN02 and 
WEN03 reconstruction using Ax = 0.05 cm and At = 0.01 s. For this experiment we constructed 
solitary waves for Boussinesq systems by solving the respective o.d.e's system in the spirit of [5] 
such that they fit to experimentally generated solitary waves before the collision. The speeds of 
the right and left-traveling solitary waves are c^^s = 0.854 m/s and ci^a = 0.752 m/s respectively. 

We observe that Boussinesq models converge to the same numerical solution with all numerical 
schemes we tested. A very good agreement with the experimental data is observed. The discrete 
mass for the Bona-Smith system is = 0.0059904310418 and for the BBM-BMM system is 
/g = 0.0059199389479 for all fluxes and reconstructions used. The variances in are mainly due 
to different types of reconstruction and not to the choice of numerical fluxes. In Table 2 these 
values are reported. 



4. Conclusions 

Initially, the finite volume method was proposed by S. Godunov [11] to compute approximate 
solutions to hyperbolic conservation laws. In the present study we made a further attempt to 
generalize this method to the framework of dispersive PDEs. This type of equations arises naturally 
in many physical problems. In the water wave theory dispersive equations have been well known 
since the pioneering work of J. Boussinesq [6] and Korteweg-de Vries [13]. Currently, the so- 
called Boussinesq-type models become more and more popular as an operational model for coastal 
hydrodynamics and other fields of engineering. 

We extend the finite volume framework to dispersive models. We tested several choices of nu- 
merical fluxes (average, Kurganov-Tadmor, characteristic), various reconstruction methods rang- 
ing from classical (TVD2, UN02) to modern approaches (WEN03, WEN05). Various choices 
of limiters have been also tested out. Advantages of speciflc methods are discussed and some 
recommendations are outlined. 
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Figure 2. Head-on collision of two solitary waves: 
Smith {9^ = 9/11), •: experimental data of [8] 
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